On présente maintenant les données que l’on a collectées et traitées afin de les utiliser dans la tâche de prédiction.
Données de comptage routier
Les données de comptage routier sont disponibles sur le site Paris Data (“Paris Data” visité le 26.12.2021), regroupant les jeux de données de la Ville de Paris, disponibles sur le site (“Données Du Comptage Routier” visité le 26.12.2021). Elles sont collectées grâce à des boucles électromagnétiques implantées dans la chaussée sur plus de 3000 tronçons de voies. L’historique des données s’étend de 2014 à ???? au pas horaire, selon la variable t_1h. Chaque boucle électromagnétique mesure le traffic sur un arc entre un noeud amont et un noeud aval. Deux données principales sont fournies : le taux d’occupation, k, qui correspond au temps de présence de véhicules sur la boucle en pourcentage d’une heure et le débit, q qui est le nombre de véhicules ayant passé le point de comptage pendant une heure. Nous gardons également d’autres variables : l’identifiant de l’arc, iu_ac, le libellé de la voie correspondante, libelle et l’état de l’arc, state. Ce dernier vaut 0 si l’état est inconnu, 1 si l’arc est ouvert à la circulation, 2 l’arc est fermé à la circulation et 3 si l’arc est invalide. Il n’est malheureusement pas expliqué la différence entre arc fermé et arc invalide dans la documentation.
Pour une année, il y a environ 29 millions observations réparties sur plus de 3000 points. On a donc procédé à une transformation des données afin de réduire leur nombre et les structurer. La première étape est l’agrégation des points d’observations selon les libellés libelle associés dans le jeu de données et à l’aide de leur identifiant iu_ac. Ensuite, on a déterminé les principaux axes de Paris en moyennant le nombre de voitures par heure des points d’observations partageant un même libellé puis en choisissant les 200 premières valeurs (hors périphérique). On obtient le graphique 1.1 et on en déduit le graphe simplifié 1.2 de Paris composé de 69 arêtes, correspondant à 69 jeux de données.
La ville de Paris rend publiquement accessible de vastes données sur la circulation routière dans la ville. La découverte du site web Paris Data (“Paris Data” visité le 26.12.2021) nous a motivé pour utiliser cette ressource dans notre projet de prédiction. Ce qui nous a surtout intrigué est la possibilité de relier les différents capteurs dans une structure de graphe. D’un côté, cela motive des études globales du réseau entier : Est-ce qu’on trouve un algorithme global, qui apprend les interactions dans le réseau entier ? De l’autre côté, l’examination de rues individuelles mène à des questions sur la corrélation locale de la circulation. Un grand boulevard se comporte-t-il comme ses avenues voisines ? Enfin, dans tout cela, nous nous poserons la question à quel horizon temporel nous réussirons à prédire la présence de voitures dans les rues de Paris.
Avant de nous lancer, nous présentons d’abord le jeu de données en détail, sa richesse mais aussi les difficultés qu’il apporte : Sa taille et les valeurs manquantes. Ayant passé une grande partie du projet sur cet aspect, nous présentons nos stratégies pour résumer et compléter ce data set.
Cela fait, nous formulons plusieurs problèmes de prédiction, à différents horizons temporels et en utilisons différentes parties du data set. (TODO: spécifier à la fin pour avoir une petite table de matières ici)
Sur le site Paris Data, on trouve énormité de jeux de données issus des activités de la ville, notamment sur le de comptage routier (“Données Du Comptage Routier” visité le 26.12.2021). Ces données sont collectées en continu grâce à des boucles électromagnétiques implantées dans la chaussée sur plus de 3000 tronçons de voies. Ces boucles et bornes de comptage mesurent principalement deux choses: Le nombre de voitures qui passent et l’occupation du tronçon de route. Disponibles de 2014 jusqu’aujourd’hui, ces deux chiffres sont mesurés au pas horaire, résultant en environ \[ 8 \text{ ans } \times 365 \text{ jours } \times 24 \text{ heures } \times 3000 \text{ capteurs } = 2.1\times10^8 \text{ lignes dans le data set.} \]
Chacune de ces lignes contient les informations suivantes qui serons plus tard nos variables explicatives \(X_i\). Les noms des variables que nous utilisons dans notre analyse plus tard sont toujours notés gras tout au long de cet exposé.
Les données historiques de circulation sont enregistrés sur le site web sous le format .txt. Chaque fichier couvre toutes les bornes pour une semaine. Pour commencer, il faut donc tout télécharger et convertir de manière automatisée ces fichiers en dataframes. En ce faisans, nous effectuons également un triage par identifiant de borne iu_ac.
En vue de la quantité énorme de 200 millions d’observations, nous sommes obligés à réduire la taille du dataset. La première démarche consiste à retirer toutes les variables inintéressantes, laissant t_1h, nbCar, rateCar et iu_ac.
Ayant réduit la dimension de chaque point dans les données, le deuxième point qu’on peut attaquer est le nombre gigantesque de senseurs de comptage. Pour ce faire, nous utilisons un modèle simplifié des rues de Paris. Comme le dataset couvre la plupart des grandes axes routières, y inclus le périphérique, notre objectif était d’agréger les données au long de ces rues :
Pour identifier où passent le plus de voitures, nous calculons des moyennes sur nbCar pour chaque libellé, c’est à dire la moyenne à travers le temps et toutes les bornes associées à chaque libellé. En prenant les 200 libellés les plus utilisées en termes de nombre de voitures (en excluant le périphérique), nous obtenons le graphique 2.1 qui représente très bien les axes routières auxquelles on s’attendrait.
Ensuite, nous ignorons les noms des libellés et regardons exclusivement ce graphique. De manière arbitraire nous en déduisons le graphe dans figure 2.2 qui représente alors un schéma très simplifié de la circulation à Paris. Grâce à notre démarche d’utiliser les rues les plus fréquentés, cette abstraction nous permet d’avoir un modèle représentatif de la circulation parisienne.
Figure 2.1: Représentation en rouge des bornes d’observation associés aux 200 libellés les plus fréquentés (hors périphérique)
Figure 2.2: Graphe simplifié. Les points noirs indiquent les intersections entre nos arêtes.
Chaque arête dans le graphe 2.2 regroupe plusieurs senseurs dont nous agrégeons les données : Après une collecte minutieuse et laborieuse des identifiants de toutes les bornes sur chacune des arêtes, nous associons pour chaque heure la moyenne de nbCar et rateCar des bornes à l’arête dont ils font partie. Cela a surtout trois avantages:
Ayant suffisamment réduit les données pour qu’ils rentrent dans le RAM d’un PC, il reste pourtant un problème: La présence de données manquantes.
En examinant des plots basiques des taux d’occupation, nous voyons souvent une image comme dans le graphique TODO. Malgré l’agrégation au long des arêtes, il reste un nombre considérable de NA dans le dataset. Afin d’executer des tâches de prédiction, il nous faut pourtant des données complètes et une partie importante de notre travail était de remplir les trous.
Figure 2.3: Taux d’occupation pour une journée ouvrière avec une valeur manquante
Figure 2.4: Visualisation de valeurs manquantes pour 30 arêtes (gris = NA)
En regardant le graphique 2.4 qu’il y a deux types de “trous” dans le data set:
Le fait que les arêtes ont de différentes structures de valeurs manquantes peut aussi être observé en regardant la répartition du pourcentage de NA’s à travers les dataframes dans figure 2.5.
Figure 2.5: Distribution des valeurs manquantes
Nous observons que les données de taux d’occupation sont généralement plus complètes : Ceci est dû au fait qu’il y a plus de capteurs de ce type. Par ailleurs, nous avons vérifié qu’une valeur manquante de nbCar implique toujours un NA dans rateCar.
Pour remplir les trous, la première étape consiste à rajouter des variables explicatives supplémentaires. En fonction du type de trou, nous verrons qu’ils auront un impact important sur la complétion.
A l’instant, nous disposons des variables nbCar et rateCar pour chaque arête et chaque heure entre 2014 et 2020. Nous exploitons deux propriétés structurelles des données pour obtenir de l’information supplémentaire:
D’abord, il s’agit de séries temporelles, donc nous rajoutons valeurs historiques de nbCar et rateCar décalés d’une heure, d’un jour et d’une semaine (nbCarLaggedHour, nbCarLaggedDay, et nbCarLaggedWeek) comme variables explicatives. L’espoir étant que la circulation se ressemble à travers des semaines consécutives par exemple. Avec le pure objectif de compléter les données manquantes, il est aussi pertinent de rajouter des valeurs du futur (rateCarFuturHour etc.), vu que le modèle d’imputation a pour but d’interpoler au lieu d’extrapoler. On a l’hypothèse que toutes ces variables décalées serviront surtout à remplir des trous locaux : Une heure ou une journée qui manque devrait facilement se laisser inférer en utilisant des données qui sont proches temporairement.
Pour les grands trous, nous n’avons souvent pas accès aux heures et jours d’avant. Cependant, grâce à l’interconnexion de nos données à travers le graphe présenté en haut, les mesures de circulation d’une arête peuvent bien être expliquées en fonction de celle de ces voisines. Dans un effort manuel, nous avons pour chaque arête \(A\) collectionné des voisins \(V_i\), dans le sens que les voitures dans \(A\) passent ensuite à l’un des \(V_i\) ou bien l’inverse. Nous faisons ici encore des choix arbitraires qui est voisin de qui, pour deux raisons. Premièrement, inférer statistiquement quelle arête influence quelle autre nécessiterait des données déjà complètes pour une régression par exemple. Deuxièmement, utiliser toutes les autres 68 arêtes au lieu de juste 2-6 voisins augmenterait trop le temps d’exécution de la complétion.
Donnons un exemple illustrant le choix de voisins ainsi que la dénomination des nouvelles variables : Parmi les voitures qui entrent dans Paris sur l’arête “pont amont - pont austerlitz” il y a une partie considérable qui continuent tout droit sur l’arête “pont austerlitz - chatelet.” Par conséquent, nous rajoutons la variable rateCar_pontausterlitzTOchatelet au dataframe de “pont amont - pont austerlitz.”
Jusqu’ici, nous n’avons utilisé que des informations déjà présentes dans le jeu de données. Dans la suite, nous rajoutons des variables extérieures qui pourraient également expliquer la circulation routière.
Variables temporelles
Le traffic routier étant relié à l’activité humaine, nous avons ajouté de nombreuses variables temporelles (notamment grâce à la librairie lubridate).
Index de la situation sanitaire en rapport avec le Covid-19
Nous avons également ajouté, covidIndex, un index allant de 0 à 100 représentant les restrictions du gouvernement sur la situation sanitaire en rapport avec le Covid-19. Il est calculé à partir de nombreux indicateurs et est fourni par l’université d’Oxford (“Covid Index” visité le 26.12.2021).
=======Index de la situation sanitaire en rapport avec le Covid-19
Nous avons également ajouté, covidIndex, un index allant de 0 à 100 représentant les restrictions du gouvernement sur la situation sanitaire en rapport avec le Covid-19. Il est calculé à partir de nombreux indicateurs et est fourni par l’université d’Oxford (“Covid Index” visité le 26.12.2021). Cette variable permettra éventuellement une étude de l’année 2020 qui voit une circulation perturbé. Pour cela, il faudra utiliser des modèles qui s’adaptent très vite à l’influence de nouvelles variables.
>>>>>>> 904223b346776040af41b6652bbfc1e13e993a53Météo
A partir de données de l’Organisation Météorologique Mondiale (“Observations Météorologiques Historiques En France” visité le 26.12.2021), nous avons extrait 2 variables météorologiques : temperature la température en Kelvin et precipitation les précipitations dans les 3 dernières heures en mm. Ces données ont été mesurées à Athis-Mons en Essonne.
On dispose d’un relevé tous les 3 heures environ donc on procède à une interpolation pour compléter les données. Etant donné que les mesures sont uniformément réparties, on utilise une interpolation linéaire basique à l’aide de la fonction na.approx de la librairie zoo.
Avec un dataframe enrichi à notre disposition pour chacune des arêtes du graphe de circulation, nous pouvons passer au remplissage des valeurs manquantes. Nous sommes pourtant restreint dans notre recherche d’un algorithme d’imputation car plusieurs variables explicatives contiennent également des NA: D’un côté, nous voulons remplir deux variables en même temps (nbCar et rateCar). De l’autre côté, les “voisins” n’ont pas moins de valeurs manquantes, laissant de grands trous dans le dataframe.
Heureusement, il y a un modèle qui peut être entraîné en dépit de valeurs manquantes parmi les covariates : Les forêts aléatoires constituées d’arbres CART. Une règle de décision dans un tel arbre peut juste être ignorée si la valeur de la variable nécessaire n’est pas renseignée. Le package miceRanger profite de ce fait en utilisant un algorithme dit d’imputation multiple : En commençant par la variable une variable \(V_1\), une forêt aléatoire est entraînée sur les lignes où \(V_1\) n’est pas manquant. Puis, les autres lignes sont “prédites” par la forêt. Avec le nouveau dataset déjà moins vide, le processus est recommence et ainsi de suite (Van Buuren 2018). Nous nous arêtons pourtant à la deuxième itération vu qu’il n’y a pas d’intérêt à compléter les données des voisins à plusieurs reprises.
Comme nous répétons le processus de complétion 69 fois, le temps d’exécution est de quelques heures sur mon laptop. Nous n’avons donc pas le luxe d’optimiser des hyperparamètres comme nous ferons plus tard pour la tâche de prédiction. Le choix de 100 arbres par forêt et 7 variables considérées à chaque split semble un bon compromis. TODO Guillaume: Voulons-nous mettre du code ici?
Au délà des capacités de complétion, les forêts aléatoires permettent de calculer des scores d’importance pour chaque variable explicatives. Ces scores peuvent aider à identifier lesquelles des variables ont le plus contribué à l’estimation de la cible (le score est calculé comme la réduction de variance à chaque split). Pour une arête du graphe, jussieu - saint michel, nous illustrons les variables qui ont le plus contribué à l’imputation :
Figure 2.6: Variables with highest importance scores for nbCar of the edge jussieu - saint michel
Figure 2.7: Variables with highest importance scores for nbCar of the edge jussieu - saint michel
Sur l’arête dont traitent les deux plots, on observe que rateCar et nbCar exhibent deux comportement assez différents : Pour le taux d’occupation, ce sont presque exclusivement les valeurs temporellement décalées qui contribuent à l’imputaion. Comme il y seulement 0.8% de valeurs à compléter, ceci confirme notre hypothèse que les petits trous d’une ou plusieurs heures sont au mieux prédites par les valeurs des heures juste avant ou après. Si nous regardons le nombre de voitures où la proportion de NA était plus importante, les deux variables ayant le plus réduit les variances aux split sont des nbCar issus d’arêtes voisines. Les trous importants en nbCar entre Jussieu et Saint Michel semblent être le mieux expliqué par ce qui se passe autour.
Un phénomène similaire peut s’observer dans les autres arêtes.(TODO Guillaume: Quoi d’autre écrire sur ceci ? Y a des scores intéressants dans RF (genre error Out of Bag) que j’aurais dû retenir de l’imputation ?)
Après l’imputation arête par arête, nous pouvons enfin compléter toutes les données grâce au fait que les variables reliés aux voisins ont été imputés dans les dataframes qui correspondent à ces voisins. Il est important de noter ici que nous allons dans la suite couper le jeu de données en deux (train - test) et que cette imputation sera effectuée individuellement sur chacune des deux parties. Faire ceci deux fois séparément est nécessaire pour préserver de l’indépendance, permettant plus tard d’estimer les erreurs de nos modèles.
Dans la suite, nous supposons que les données sont complètes. Il faut pourtant remarquer que ces forêts aléatoires font indirectement partie des modèles de prédiction que nous allons utiliser.
Choix des jeux de données
On découpe notre jeu de données en 2 parties de proportion 2/3 et 1/3 : une partie apprentissage de 2014 à 2017 et une partie test de 2018 à 2019. Ainsi, on pourra évaluer les performances de nos modèles sur les différentes périodes d’une année.
Choix de la métrique
Pour évaluer ces performanes, on choisit la métrique Root Mean Square Error (RMSE) qui représente la racine carrée du deuxième moment d’échantillonnage des différences entre les valeurs prédites et les valeurs observées. On n’utilise pas la métrique Mean Absolute Percentage Error (MAPE) car un nombre important des valeurs de rateCar sont nulles.
Modèles naïfs à battre
Afin de comparer nos modèles, on construit 3 modèles témoins. Ces modèles sont dit naïfs car extrêmement simpliste :
naiveModel1 prévoit le taux d’occupation d’une heure \(t \in \{ 0, \dots, 23 \}\) d’un jour \(j \in \{ 1, \dots, 7\}\) d’un mois \(m \in \{ 1, \dots, 12 \}\) d’une année \(a \in \{ 2018, 2019 \}\) en moyennant les taux d’occupation de l’heure \(t\) du jour \(j\) du mois \(m\) pour toutes les années \(a \in \{ 2014, \dots, 2017 \}\).
naiveModel2 fait de même en calculant la médiane et non la moyenne.
naiveModel3 est simplement le taux d’occupation à l’heure précédente rateCarLaggedDay.
On calcule la moyenne des erreurs RMSE pour les 69 arêtes :
| naiveModel1 | naiveModel2 | naiveModel3 |
|---|---|---|
| 5.844 | 5.899 | 3.955 |
Recherche de paramètres
<<<<<<< HEADAfin d’optimiser les performances de nos différents modèles, nous effectuons de la recherche de paramètres. On utilise les 2 méthodes suivantes :
Validation croisée sur 16 blocs. On divise ainsi nos 4 années de données d’apprentissage en blocs de 3 mois, ce qui englobe les différentes saisonnalités de nos données. On a codé cette méthode à la main dans la fonction ….
Validation croisée progressive, où on fixe une fenêtre initiale d’apprentissage de 2 ans, que l’on incrémente dans l’ordre chronologique de 2 mois à chaque itération, pour entraîner le modèle et mesurer sa performance à l’horizon 1 (c’est-à-dire prédire le taux d’occupation de l’heure suivante). Pour utiliser cette méthode, on utilise le package caret avec le paramètre timeSlice.
Figure 1.4: Validation croisée progressive (“Validation Croisée Progressive” visité le 03.02.2022)
La caractéristique des données de type série temporelle est qu’elles sont ordonnées et possèdent une dépendance temporelle. C’est pourquoi il semble pertinent de ne pas mélanger les observations avant de découper le jeu de données en données apprentissage et données test. De plus, il semble aussi important de ne pas utiliser des données futures pour prédire le passé. Dans la validation croisée, il est donc naturel de retirer les blocs d’apprentissage ultérieurs au bloc test. C’est la validation croisée progressive. Cependant, en pratique, on souhaite entraîner le modèle une unique fois sans le mettre à jour tous les 2 mois. Dans ce cas, la validation croisée est plus pertinente.
=======Afin d’optimiser les performances de nos différents modèles, nous effectuons de la recherche de paramètres à l’aide du package caret. On utilise les 2 méthodes suivantes :
cv : validation croisée sur 16 blocs. On divise ainsi nos 4 années de données d’apprentissage en blocs de 3 mois, ce qui englobe les différentes saisonnalités de nos données.
timeSlice : où on fixe une fenêtre initiale d’apprentissage de 2 ans, que l’on incrémente dans l’ordre chronologique de 2 mois à chaque itération, pour entraîner le modèle et mesurer sa performance à l’horizon 1 (c’est-à-dire prédire le taux d’occupation de l’heure suivante).
L’avantage de timeSlice est qu’il permet d’optimiser un modèle pour la prédiction à un horizon déterminé. De plus, contrairement à la validation croisée, il n’utilise pas de données futures pour prédire le passé. Cependant, en pratique, on souhaite entraîner le modèle une unique fois sans le mettre à jour tous les 2 mois. Dans ce cas, la validation croisée est plus pertinente.
>>>>>>> 904223b346776040af41b6652bbfc1e13e993a53Lors de nos recherches de paramètres, nous avons observé que les 2 méthodes sélectionnent des paramètres proches. Finalement, on se base sur la validation croisée en 16 blocs pour la recherche de paramètres.
Pour permettre la reproductibilité des résultats, on fixe la graine du générateur aléatoire. De plus, on effectue la recherche de paramètre pour une seule arête sur les 69 car c’est un processus long et on fait l’hypothèse que le paramètre optimal pour la prévision sur une arête est proche de celui de toutes les autres arêtes.
Figure 2.1: Recherche du paramètre cp optimal =======
Figure 3.1: Recherche du paramètre cp optimal >>>>>>> 904223b346776040af41b6652bbfc1e13e993a53